function vns(inDir,outDir)

%inDir='C:\Users\matthew\Desktop\workspace\vnsoddmc\data\vns\';
%outDir='C:\Users\matthew\Desktop\workspace\vnsoddmc\analysis\vns\';

filter=fullfile(inDir,'*_*');
subjs=dir(filter);
nSubjs=length(subjs);
nSamples=1024;

Tr2use=minTrial2use(inDir,nSamples);


for i=1:nSubjs
    filename=fullfile(inDir,subjs(i).name);
    data=importEDF2FT(filename);
    trials=rechunk(data,nSamples,Tr2use);
    [outfile,r]=strtok(subjs(i).name,'.');
    outfile=strcat(outfile,'-trials');
    save(fullfile(outDir,outfile),'trials')
end

pliXBand(outDir,outDir,Tr2use);
pliStat(outDir,outDir)




function pliXBand(inDir,outDir,Tr2use)

bands={[1 4],[4 8],[8 13],[13 22],[22 49]};
bandLabels={'delta','theta','alpha','beta','gamma'};
nBands=length(bandLabels);

filter=fullfile(inDir,'*-trials*');
files=dir(filter);
nSubjs=length(files);


for i=1:nSubjs
    
    subjName=files(i).name;
    load(fullfile(inDir,subjName),'trials')
    subjName=strtok(subjName,'-');
    for j=1:nBands   
        
        bf_trials=selectBand(trials,bands{j})        
        pliXT=pliXTrial(bf_trials,Tr2use);
        
        outfile=strcat('pli_',bandLabels{j});
        outfile=strcat(outfile,'-');
        outfile=strcat(outfile,subjName);
        
        save(fullfile(outDir,outfile),'pliXT');
    end    
end


function PLI=pli(a)

% a is a filtered multichannel signal (time x channels)
% hilbert(a) calculates analytic signal (complex valued) of each
% column of a. Phase Lag Index between channel i and j averaged over
% time bins is stored in PLI(i,j)
% number of channels
N=size(a,2);
PLI(1:N,1:N)=0;
complex_a=hilbert(a);
for i=1:N-1
   for j=i+1:N
        PLI(i,j)=abs(mean(sign(imag(complex_a(:,i)./complex_a(:,j)))));
   end
end
PLI=PLI+PLI';

function [pli3d]=pliXTrial(bf_trials,nTrials2use)
if(nTrials2use==-1)
    nTrials2use=length(bf_trials.trial);
end
nChs=size(bf_trials.trial{1},1);
pli3d=zeros(nTrials2use,nChs,nChs);
for i=1:nTrials2use
   pli3d(i,:,:)=pli(bf_trials.trial{i}');  
end



function pliStat(inDir,outDir)

lf=[1 3 4];
lc=9;
lt=[8 13];
lp=14;
lo=18;

rf=[2 6 7];
rc=11;
rt=[12 17];
rp=16;
ro=19;

zz=[5 10 15];

areas={lf;lc;lt;lp;lo;rf;rc;rt;rp;ro;zz};
bandlabels={'delta','theta','alpha','beta','gamma'};
alphalevel=0.05;

for i=1:length(bandlabels)
    
    filter=strcat('pli_',bandlabels{i});
    filter=strcat(filter,'*_A*');
    subjs=dir(fullfile(inDir,filter));

    nSubjs=length(subjs);
    for j=1:nSubjs

            fileName=subjs(j).name;
            load(fullfile(inDir,fileName),'pliXT');
            one_pliXT=pliXT;
        
            fileName(end-4)='D';

            load(fullfile(inDir,fileName),'pliXT');

            two_pliXT=pliXT;

            
            statsingle=statOverTr(one_pliXT,two_pliXT,areas,alphalevel);
            
            [outfile,r]=strtok(fileName,'.');
            [r,outfile]=strtok(outfile,'_');
            outfile=strcat('stat',outfile);
            outfile=outfile(1:end-1);
            outfile=strcat(outfile,'AD');
             
            save(fullfile(outDir,outfile),'statsingle');
         
    end
    
    statOverSubjs(inDir,outDir,bandlabels{i},alphalevel)
end




function [trials]=rechunk(data,nSamples,Tr2use)
cfg=[];

len=nSamples/data.fsample;
cfg.length=len;  
cfg.overlap=0;
trials=ft_redefinetrial(cfg,data);
if(Tr2use~=-1)
    trials.trial=trials.trial(1:Tr2use);
    trials.time=trials.time(1:Tr2use);
end

function [bf_trials]=selectBand(trials,band)

bf_trials=trials;
nTrials=length(trials.trial);

for i=1:nTrials
    
    bf_trials.trial{i}=eegfilt(trials.trial{i},trials.fsample,band(1),band(2));

end
function [cdata]=importEDF2FT(filename)

chLabelslOI={'MKR+','Fp1','Fp2','F7','F3','Fz','F4','F8','T3','C3','Cz','C4','T4','T5','P3','Pz','P4','T6','O1','O2'}
%filename='C:\Users\matthew\Desktop\Projects\VNS\data\Autiero 5.05.08 7y.edf';
%channel of interest
to_import=pop_biosig(filename);
chLabels={to_import.chanlocs.labels};
choi=zeros(1,length(chLabels));
for i=1:length(chLabelslOI)
    choi=choi+strcmp(chLabelslOI{i},chLabels);
end
indxChOI=1:length(chLabels);
choi=logical(choi);
data.trial{1}=to_import.data(indxChOI(choi),:);
indxMKR=strcmp('MKR+',chLabels);
choi(indxMKR)=0;%mk+
data.label={to_import.chanlocs(indxChOI(choi)).labels}';



data.label

data.fsample=to_import.srate;

cdata=chunkData(data);
function [cdata]=chunkData(data)

cdata=data;   
mk=data.trial{1}(1,:);
mk=abs(mk);
indx=mk>1;
indxv=1:length(indx);
v=indxv(indx);
indx=splitVec(v,'consecutive');

[s ntrial]=size(indx);
cdata.trial=cell(s,ntrial);
sampl_interval=1/cdata.fsample;

for i=1:ntrial
    cdata.trial{i}=data.trial{1}(2:end,indx{i});
    cdata.time{i}=(0:1:length(cdata.trial{i})-1)*sampl_interval;
end

function [minTr2use]=minTrial2use(inDir,trialLen)

filter=fullfile(inDir,'*_*');
subjs=dir(filter);
nSubjs=length(subjs);

minNTrials=zeros(1,nSubjs)
for i=1:nSubjs
    %subjsNames{i}=subjs(i).name
    filename=fullfile(inDir,subjs(i).name);
    data=importEDF2FT(filename);
    trials=rechunk(data,trialLen,-1);
    
    minNTrials(i)=length(trials.trial);
end
minTr2use=min(minNTrials);




function [res]=statOverTr(one_pliXT,two_pliXT,areas,alphalevel)

nTrial=min(size(one_pliXT,1),size(two_pliXT,1));
%one_trialXintra_area=zeros(nTrial,length(areas));
%two_trialXintra_area=zeros(nTrial,length(areas));

interAreaLen=((length(areas)*length(areas))-length(areas))/2;

one_trialXinter_area=zeros(nTrial,interAreaLen);
two_trialXinter_area=zeros(nTrial,interAreaLen);


s=size(one_pliXT);
indx=ones(s(2),s(3));
indx=indx-tril(indx);
indx=logical(indx);

one_avg_pliXtrial=zeros(nTrial,1);
two_avg_pliXtrial=zeros(nTrial,1);


for i=1:nTrial
   
   tone=squeeze(one_pliXT(i,:,:));
   ttwo= squeeze(two_pliXT(i,:,:));
    
 %  one_trialXintra_area(i,:)=shortPLI(tone,areas);
  % two_trialXintra_area(i,:)=shortPLI(ttwo,areas);
   
   one_trialXinter_area(i,:)=longPLI(tone,areas);
   two_trialXinter_area(i,:)=longPLI(ttwo,areas);
   
   one_avg_pliXtrial(i)=mean(tone(indx));
   two_avg_pliXtrial(i)=mean(ttwo(indx));

end

%res.intra={one_trialXintra_area,two_trialXintra_area};
res.inter={one_trialXinter_area,two_trialXinter_area};
res.avg={one_avg_pliXtrial,two_avg_pliXtrial};

%[res.h_intra,res.p_intra]=ttest(one_trialXintra_area,two_trialXintra_area,alphalevel);
[res.h_inter,res.p_inter]=ttest(one_trialXinter_area,two_trialXinter_area,alphalevel);
[res.h_avgpli,res.p_avgpli]=ttest(one_avg_pliXtrial,two_avg_pliXtrial,alphalevel);

function statOverSubjs(inDir,outDir,bandName,alphalevel)

filter=strcat('stat_*',bandName);
filter=strcat(filter,'*');

files=dir(fullfile(inDir,filter));
nSubjs=length(files);

load(fullfile(inDir,files(1).name),'statsingle');

%nAreas=size(statsingle.intra{1},2)
nPairAreas=size(statsingle.inter{1},2);

%one_SubjXintraA=zeros(nSubjs,nAreas);
%two_SubjXintraA=zeros(nSubjs,nAreas);
one_SubjXinterA=zeros(nSubjs,nPairAreas);
two_SubjXinterA=zeros(nSubjs,nPairAreas);
one_SubjAVG=zeros(nSubjs,1);
two_SubjAVG=zeros(nSubjs,1);

for i=1:nSubjs
    load(fullfile(inDir,files(i).name),'statsingle');

    %one_SubjXintraA(i,:)=mean(statsingle.intra{1},1);
    %two_SubjXintraA(i,:)=mean(statsingle.intra{2},1);
    one_SubjXinterA(i,:)=mean(statsingle.inter{1},1);
    two_SubjXinterA(i,:)=mean(statsingle.inter{2},1);
    one_SubjAVG(i)=mean(statsingle.avg{1},1);
    two_SubjAVG(i)=mean(statsingle.avg{2},1);
    
    [r,subjName]=strtok(files(i).name,'-');
    subjName=subjName(2:end);
    subjName=strtok(subjName,'.');
    multi.subjNames{i}=subjName;  
end

%multi.intra={one_SubjXintraA,two_SubjXintraA};
multi.inter={one_SubjXinterA,two_SubjXinterA};
multi.avg={one_SubjAVG,two_SubjAVG};

%[multi.h_intra,multi.p_intra]=ttest(one_SubjXintraA,two_SubjXintraA,alphalevel);
[multi.h_inter,multi.p_inter]=ttest(one_SubjXinterA,two_SubjXinterA,alphalevel);
[multi.h_avgPLI,multi.p_avgPLI]=ttest(one_SubjAVG,two_SubjAVG,alphalevel);

outfile=fullfile(outDir,strcat('statMultiSubjs_',bandName));

save(outfile,'multi');
% function [v]=shortPLI(pli,areas)
% 
% v=zeros(1,length(areas));
% for i=1:length(areas)
%     m=pli(areas{i},areas{i}); 
%     indx=ones(size(v));
%     indx=indx-tril(indx);
%     indx=logical(indx);    
%     v(i)=mean(m(indx));
% end
function [interArea]=longPLI(pli,areas)

v=zeros(length(areas),length(areas));
for i=1:length(areas)-1
    for j=i+1:length(areas)
        m=pli(areas{i},areas{j});
        s=size(m);
        v(i,j)=mean(reshape(m,s(1)*s(2),1));
    end
end

indx=ones(size(v));
indx=indx-tril(indx);
indx=logical(indx);
%indx=logical(indx);%%warning if phase is 0
interArea=v(indx)';






